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Abstract. The equations governing the vertical struc- 
ture of a stationary keplerian accretion disc supporting 
an Eddington atmosphere are presented. The model is 
based on the a-prescription for turbulent viscosity (two 
versions are tested) , includes the disc vertical self-gravity, 
convective transport and turbulent pressure. We use an 
accurate equation of state and wide opacity grids which 
combine the Rosseland and Planck absorption means 
through a depth-dependent weighting function. The nu- 
merical method is based on single side shooting and in- 
corporates algorithms designed for stiff initial value prob- 
lems. A few properties of the model are discussed for a 
circumstellar disc around a sun-like star and a disc feed- 
ing a 10 8 Mq central black hole. Various accretion rates 
and a-parameter values are considered. 

We show the strong sensitivity of the disc structure 
to the viscous energy deposition towards the vertical axis, 
specially when entering inside the self-gravitating part of 
the disc. The local version of the a-prescription leads to a 
"singular" behavior which is also predicted by the verti- 
cally averaged model: there is an extremely violent density 
and surface density runaway, a rapid disc collapse and a 
temperature plateau. With respect, a much softer transi- 
tion is observed with the "a'P-formalism" . Turbulent pres- 
sure is important only for a > 0.1. It lowers vertical den- 
sity gradients, significantly thickens the disc (increases its 
flaring) , tends to wash out density inversions occurring in 
the upper layers and pushes the self-gravitating region to 
slightly larger radii. Curves localizing the inner edge of the 
self-gravitating disc as functions of the viscosity parame- 
ter and accretion rate are given. The lower a, the closer to 
the center the self-gravitating regime, and the sensitivity 
to the accretion rate is generally weak, except for a > 0.1. 

This study suggests that models aiming to describe 
T-Tauri discs beyond about a few to a few tens astronom- 
ical units (depending on the viscosity parameter) from the 
central protostar using the cv-theory should consider ver- 
tical self-gravity, but additional heating mechanisms are 
necessary to account for large discs. The Primitive Solar 
Nebula was probably a bit (if not strongly) self-gravitating 



Send offprint requests to: Jean-Marc.Hure@obspm.fr 



at the actual orbit of giant planets. In agreement with 
vertically averaged computations, a-discs hosted by ac- 
tive galaxies are self-gravitating beyond about a thousand 
Schwarzchild radii. The inferred surface density remains 
too high to lower the accretion time scale as requested to 
fuel steadily active nuclei for a few hundred millions years. 
More efficient mechanisms driving accretion are required. 

Key words: Accretion, accretion disks - Equation of 
state - Solar system: formation - Galaxies: active - Galax- 
ies: nuclei 



1. Introduction 

Viscous discs can acquire large dimensions under the effect 
of angular momentum redistribution. This is corroborated 
by the observation of Young Stellar Objects (YSO) which 
reveals the emission of wide gaseous discs with an outer 
radius reaching 100 - 1000 AU (Beckwith et al., 1990; 
Pudritz et al., 1996; Duvert et al., 1998; Guilloteau & 
Dutrey, 1998; Shepherd & Kurtz, 1999). Also, discs hosted 
by Active Galactic Nuclei (AGN) , although yet unresolved 
— except in the active galaxy NGC4258 (Herrnstein et al., 
1999) — could probably be of much larger size. The exis- 
tence of strong similarities between YSOs and AGNs, like 
jets and outflows (Falcke, 1998; Frank, 1998), some ex- 
cesses in the spectral energy distribution at infrared wave- 
lengths (Adams & Shu, 1986; Zdziarski, 1986; Bertout, 
1989; Sanders et al., 1989; Voit 1991; Kenyon, Yi & Hart- 
mann, 1996), a surrounding (dusty) torus (Gusten, Chini 
& Neckel, 1984; The k Molster, 1994; Drinkwater, Combes 
& Wiklind, 1996; Sandqvist, 1999) indicates that the same 
physical mechanisms should take place in discs, despite 
the difference in the central mass scale. This constitutes 
an interesting challenge for the disc theory, especially to 
understand the transport of angular momentum. 

Many models are based on the a-theory of thin discs 
(Shakura & Sunyaev, 1973; Pringle, 1981) and assume 
a vertically averaged structure (e.g., Collin & Dumont 
1990; Ruden & Pollack, 1991; Cannizzo & Reiff, 1992; 
Hure et al. 1994a; Artemova et al. 1996; Burderi, King 
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& Szuszkiewicz, 1998; Drouart et al. 1999; Aikawa et al., 
1999). A physically more satisfactory approach is the in- 
vestigation of the vertical structure in detail, as done 
for stars. Such a problem has been discussed by several 
groups already, in different contexts, with various goals 
and degrees of sophistication: at the scale of AGNs (Can- 
nizzo, 1992; Wehrse, Storzer & Shaviv, 1993; Dorrer et al., 
1996; Siemiginowska, Czerny & Kostyunin, 1996; Hubeny 
& Hubeny, 1998 and references therein; Sincell & Kro- 
lik, 1997; Rozanska et al., 1999), in Cataclysmic Variables 
(CVs) (Smak, 1984; Mineshige & Osaki, 1983; Meyer & 
Meyer-Hofmeister, 1982; Pojmahski, 1986; see Cannizzo, 
1993 for a review; Milsom, Chen & Taam; 1994; Dubus et 
al., 1999) and in YSOs (Malbet & Bertout, 1991; Bell et 
al., 1997; D'Alessio et al., 1998) including the Primitive 
Solar Nebula (PSN) (Lin & Papaloizou, 1980; Papaloizou 
& Terquem, 1999). Basic considerations show that the 
outermost regions of (low mass) discs are expected to be 
regulated by self-gravity (Goldreich & Lynden-Bell, 1965; 
Shlosman & Begelman, 1987). This may be the case of 
AGN and YSO discs. A few constraints on the disc thick- 
ness, mass and rotation motion at large radii are set by 
observations (Guilloteau & Dutrey, 1998; Mundy, Looney 
& Welch, 2000; Herrnstein et al., 1999). So, it appears es- 
sential to construct models for as accurate and realistic as 
possible, despite the lack of knowledge regarding turbu- 
lent viscosity which causes accretion. To our knowledge, 
no self-consistent 2D-model accounting for vertical self- 
gravity in the outermost parts of discs has been published 
yet. This is the aim of this paper. The present model in- 
cludes simultaneously vertical convection, self-gravitation, 
turbulent pressure, and realistic equation of state (EOS) 
and opacities. Of particular interest here is the transition 
from the "classical" disc to the self-gravitating disc which 
is predicted to be as close as a few 10 3 — 10 4 i?* (i?* is the 
radius of the central object) by vertically averaged models 
(Ruden & Pollack, 1991; D'alessio, Calvet & Hartmann, 
1997; Hure, 1998). The hypothesis of the model and rele- 
vant equations are developed in Sect. 2. Two versions of 
the a-prescription are tested. The ingredients (EOS and 
opacities) as well as the numerical method are presented 
in Sect. 3. We discuss in Sect. 4 a few properties of the 
model, namely the effect of turbulent pressure and depth 
dependent viscosity, and the position of the inner edge of 
the self-gravitating disc, for prototypal circumstellar and 
AGN discs, for various accretion rates. An Appendix con- 
tains a note concerning the treatment of convection and a 
high precision formula fitting the EOS. 

2. Model for the vertical structure: hypothesis 
and relevant equations 

2.1. General considerations. Accounting for self-gravity 

A distinctive feature of steady state keplerian accretion 
discs is the absence of coupling between the vertical struc- 
ture and the radial structure (e.g. Frank, King & Raine, 



1992). This attractive property is undoubtedly an over- 
simplification and probably does not match reality. How- 
ever, it is particularly advantageous from a modeling point 
of view because any annulus can be treated individually, 
whatever the state of neighbors, unlike in thick discs where 
pressure gradients and energy transport in the radial di- 
rection gain in importance with respect to their vertical 
counterparts (Maraschi, Reina & Treves, 1976; Abramow- 
icz, Calvani & Nobili, 1980; Abramowicz et al., 1988; 
Narayan, Madevan & Quataert, 1998). Note that the kep- 
lerian assumption which fixes the rotation law to the value 

where M is the central mass and R the polar radius, re- 
quires that the gas remains confined at altitudes z such as 
z 2 <C R 2 , and small radial pressure gradients too. 

Self-gravity may influence and even dominate the equi- 
librium structure and dynamical evolution of almost any 
kind of disc, not only in massive or thick discs or tori 
where effects are global (Bodo & Curir, 1992; Hashimoto, 
Eriguchi, Muller, 1995; Boss, 1996; Laughlin & Rozyczka, 
1996, Masuda, Nishida & Eriguchi, 1998), but also in low 
mass keplerian discs (Paczynski, 1978; Shlosman & Begel- 
man, 1987) as soon as the mass density of the accreted gas 
exceeds ~ f2 2 /AnG locally, f2 being the rotation frequency. 
In that latter case which is of interest here, self-gravity 
increases vertical pressure gradients and gathers matter 
closer to the midplane (Sakimoto & Coroniti, 1981; Shore 
& White, 1982; Cannizzo & Reiff, 1992; Hure et al, 1994; 
Hure, 1998). So, only the outermost regions of keplerian 
discs where f2 reaches low values can be affected by vertical 
self-gravity, except in very special situations (e.g. Sincell 
& Krolik, 1997). 

Accounting correctly for the disc gravity requires the 
resolution of the Poisson equation (Hunter, 1963; Storzer, 
1993; Bertin & Lodato, 1999) 

A$ disc (i?, z) = -AttG P {R, z) (2) 

where p is the gas mass density and $ dlsc is the gravita- 
tional potential due to the bare disc (the </>-invariance is 
assumed). This is a very difficult task, in particular be- 
cause the disc surface has a non trivial form which is not 
known a priori and the deviation from sphericity is ex- 
treme. As Eq.(|^) connects all annuli together, the method 
of "independent rings" no longer applies, unless an itera- 
tive scheme in which the potential would be step by step 
improved from the actual density field until convergence 
(Stahler, 1983). There is however no guarantee that such 
a scheme effectively converge and the computational time 
might be prohibitive (Eriguchi & Muller, 1991). Here, we 
follow another, more simple approach: we adopt the infi- 
nite and i?-homogeneous slab approximation (Paczynski, 
1978; Sakimoto & Coroniti, 1981; Shore & White, 1982; 
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Cannizzo & Reiff, 1992; Liu, Xie & Ji, 1994; Hure, 1998) 
which yields the gravity due to the disc 



where A„ = 



<9$ 



disc 



dz 



-47rGS(i?, z), 



where the surface density E is defined by 



E(i?,z) 



p(R, z')dz' , 



(3) 



(4) 



but other assumptions are possible (Mineshige & 
Umemura, 1997). This enables again to investigate the 
disc structure annulus by annulus, but introduces a bias: 
it tends to overestimate (underestimate) self-gravity in 
regions of high (respectively low) densities. We expect 
that the present approximation gives acceptable results, 
at least in regions where the surface density radial gradi- 
ents remain low with respect to density. Anyway, the error 
made with respect to a self-consistent model is not known 
and would be interesting to estimate. 

2. 2. The equations for the disc interior 

When the upward transport of heat is treated in the dif- 
fusion approximation, the vertical structure of a keplerian 
disc, possibly irradiated, can be determined from the res- 
olution of a system of four first order coupled ordinary 
differential equations (ODEs) in between the midplanc 
and the top of the disc (Pojmanski, 1986; Tuchman, Mi- 
neshige & Wheeler, 1990; Cannizzo, 1992; Meyer & Meyer- 
Hofmeister, 1982; Milsom, Chen & Taam, 1994; Dubus 
et al., 1999). The complexity of the problem rises when 
a multi-frequency radiative transfer is performed, as re- 
quired to compare theoretical spectra with observations 
and make key predictions (Ross, Fabian & Mineshige, 
1992; Wehrse, Storzer & Shaviv, 1993; Dorrer et al., 1996; 
Sincell & Krolik, 1997; D'alessio et al. 1998; El-Khoury & 
Wickramasinghe, 1998; Hubeny & Hubeny, 1998; De Kool 
& Wickramasinghe, 1999). The four equations specify re- 
spectively (Frank, King & Raine, 1992): 

— the pressure gradient V Z P which describes the hydro- 
static equilibrium of each fictitious slab 



1 dP 

p dz 



= -n 2 z-4irGY; = g z 



(5) 



— the heat flux gradient V Z F due to viscous heating 

dF 9 a 
— = -pvil , 

dz 4 



(G) 



where v is the z-dependent viscosity law (see Sect. 2.4), 
the temperature gradient V Z T which is determined by 
the heat net flux transported upwards through radia- 
tion and convection 



dT 

dz 



-T 



(J) 



jj^p is the pressure height scale and 
V = jp-p is the actual gradient (see the Appendix, 
Sect. A), 

— the surface density gradient V Z E 



dE _ 
dz ^ 



(8) 



Note that this last equation is not relevant when self- 
gravity is left apart (the total surface density of the disc 
+ atmospheres E t = 2 x E(i?,oo) can be easily computed 
a posteriori from Eq. (|4j) , once the mass density distribu- 
tion is known). The above ODEs must be supplemented 
by a closure relation, an equation of state. For a mixture 
of radiation and perfect gas undergoing atomic ionization 
and molecular dissociation at LTE, the total pressure is 
linked to the density and temperature of matter by 



P 



pkT A a 
pmn 3 c 



T 



(9) 



where /imn is the mean mass per particle which depends 
on the density and temperature (see Sect. 3.1 and the 
Appendix, Sect. B). The above expression for radiation 
pressure is compatible with an optically thick disc only. 

2.3. Accounting for turbulent pressure in the framework 
of the a-prescription 

Turbulence is the main mechanism driving accretion in 
discs. It is an extra source of pressure. According to the 
standard theory of a-discs (Shakura & Sunyaev, 1973; 
Pringle, 1981), the typical velocity of turbulent eddies can 
be taken as ~ \/ac s if one assumes the equipartition be- 
tween length and velocity turbulent scales (c s is the adia- 
batic sound speed). So, the turbulent pressure is 



p t = aTiP 



(10) 



where T\ = ( ^-j~ )ad is the first adiabatic exponent in the 
total (gas plus radiation) pressure. The effect of turbulent 
pressure on the hydrostatic equilibrium is therefore ex- 
pected if a is not too small, as simulations shall confirm. 
If we include pt into Eq.(|^) and rewrite it in terms of a 
density gradient equation, we find 



dz P K > Xp (1 



(11) 



where \t and \p are respectively the temperature and 
density exponents of the total pressure, and we have as- 
sumed VJi = (this is justified since 0.9 < Y\ < 1.7, 



see for instance the Appendix, Fig.(A.l)). This expression 
clearly shows the existence of a density inversion (that 
is, a zone where \7 z p > 0) each time the actual gradient 
satisfies xtV > 1 (Rozanska et al., 1999). Such an in- 
version is likely to occur in radiative pressure dominated 
layers where xt is the largest or/and in convectively un- 
stable zones where V may be large. Note that the inversion 
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still remains if turbulent pressure plays a role, but with a 
weaker amplitude. 

It is likely that turbulent pressure plays a role, not 
only on the hydrostatic equilibrium as considered here, 
but also on advection of matter and energy both radially 
and vertically. We ignore these effects. 

2.4- On depth- dependent viscosity laws 

In vertically averaged disc models, the anomalous vis- 
cosity takes locally a mean value defined from midplane 



quantities, namely v 



— ^midplanef 



(= ac s 2 (0)/ft). Here, we 



need to specify how v varies with the altitude. This point 
is critical since our knowledge of turbulent viscosity re- 
mains extremely limited, if not null. There are however 
two standard ways to do this within the framework of the 
a-prescription. The most common one is to assume that 
the shear stress is proportional to some pressure (the so- 
called " a'P-formalism" ; e.g. Cannizzo, 1992). It leads to a 
viscosity law 



2aP 
3 ftp 



Vi(R, z) 



(12) 



where V can be either gas pressure, or total pressure 
or some combination of the two (Siemiginowska, Czerny 
& Kostyunin, 1996; Artemova et al., 1996; Hameury et 
al., 1998; Papaloizou & Terquem, 1999), with different 
consequences on the disc stability (Camenzind, Demolc 
& Straumann 1986; Clarke, 1988). The other way uses 
the local version of the a-prescription (Meyer & Meyer- 
Hofmeister, 1982) 



v 2 



ac s \p 



v 2 (R,z) 



(13) 



where X p is the reduced pressure height scale (see the Ap- 
pendix, Sect. A). Whatever the expression the authors 
adopt, the a-parameter is most often assigned to a fixed 
value in a disc. However, theoretical arguments, numerical 
simulations as well as observational constraints indicate 
that a should vary, not only with the radius (Shakura & 
Sunyaev, 1973; Cannizzo, 1993; Lasota& Hameury, 1998), 
but also with the altitude via the temperature or density, 
or other quantities (Brandenburg, 1998). Strictly speak- 
ing, the need for strong variations of the a-parameter 
means the failure of the a- viscosity model. 

Since v\ and v 2 are formally different, even with the 
same a-parameter, they lead to different vertical struc- 
tures and consequently to different discs, as we shall see 
below. Besides, the correspondence between the two laws 
(a multiplying factor |ftA p ri/c s from V\ to v 2 , if V = P) 
is much more subtile than shifting a: in general, ftA p =/= c s 
at any altitude (this is specially true at the equatorial 
plane where A p — > oo; see Sect. 4.3). The volumic energy 
production associated to v 2 is 



dz 4 * 



(14) 



Note that both Eqs.@ and @ satisfy || < above the 
equatorial plane, meaning that the gas is accreted much 
faster at the equatorial plane than at the top of the disc. 
So, these laws implicitly suggest the existence a plan par- 
allel shear which might be able to trigger turbulence and 
to expand turbulent eddies in the radial direction. 

There are still no physical arguments to decide if v 2 is 
better than v\ . That is why we use both expressions in the 
following (in particular, for v\ we take V = P). In fact, a 
wide class of functions v(R, z) should be considered and 
their effects compared. It is possible that turbulent viscos- 
ity shows a weaker dependence with z than considered so 
far and even does not depend on the altitude at all. This is 
precisely the case with the (5- viscosity prescription which 
is suggested by laboratory experiments (Pringle & Rees, 
1972; Lynden-Bell & Pringle, 1974; Richard & Zahn, 1999; 
Duschl, Bicrmann & Strittmatter, 2000; Hure, Richard & 
Zahn, 2000). 

2. 5. Equations for the Eddington atmosphere 

In the Eddington approximation, the structure of the at- 
mosphere is governed by three (only two in the absence 
of self-gravity) first orders ODEs. Two of these (namely 
Eqs.(||) and (g)) are basically the same as for the disc in- 
terior, and the third one controls the optical depth r in 
the atmosphere 



dT 

dz 



(15) 



where k is a grey absorption coefficient. In order to 
prevents any temperature runaway which would lead to 
the formation of a hot corona (Shaviv & Wehrse, 1991; 
Rozanska et al., 1999), we assume that there is no viscous 
energy generation in this layer and no active turbulence. It 
means that the atmosphere is not accreted. We are aware 
that the validity of all these approximations, including the 
Eddington approximation, is easily open to criticism. The 
temperature within the atmosphere follows the law (Mi- 
halas, 1978) 



rp4 



off 



T+ 3 



(16) 



where r < | and T e g is the effective temperature which 
is fixed by the outgoing flux due to internal and external 
heating sources, at the photosphere's base. Actually, the 
global effect of the disc irradiation can be taken into ac- 
count at this level by a suitable definition of T g (Kenyon 
& Hartmann, 1987; Hubeny, 1990; Robinson, Marsh & 
Smak, 1993; Dubus et al, 1999). Wee see from Eqs.(§), 
( |l5| ) and ( |l6| ) that the density gradient is 



dp 

dz 



(ft 2 z + 47rGE) 



XT up 



4 Xp (r- 



(17) 



and may also become positive. As suggested above, we do 
not take into account turbulent pressure in this layer. This 
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induces a discontinuity in the density (or gas pressure) 
gradient at the altitude where the disc and its atmosphere 
join together, but not in the density itself. 

2.6. Boundary conditions and matching conditions at the 
disc/ atmosphere interface 

There are two problems to solve simultaneously: the struc- 
ture of the disc interior from Eqs.(^)(^)(||) and (p"l| ) and 
the structure of its atmosphere from Eqs.(P)(fL5|) and 
(p"7j). These differential equations are subject to Dirich- 
let boundary conditions. At the midplane, the symmetry 
imposes S = and F = like for a star. At the top 
atmosphere located at z = H (to be determined), con- 
ditions are P = P am b the ambient pressure (which may 
include the pressure due to external illumination), r = 
and £ = I St, half the total surface density of the disc 
and atmosphere. Some authors specify the density (for in- 
stance p — ► 0) instead of the pressure at the boundary (e.g. 
Pojmahski, 1986; Dorrer et al. 1993; El-Khoury & Wick- 
ramasinghe, 1999). We choose a value typical of the inter- 
stellar medium, namely P am b/k — 10 5 K.cm~ 3 (Duley & 
Williams, 1988), significantly smaller than in D'alessio et 
al. (1998). At the base of the atmosphere located at z = h 
(to be determined too), the matching conditions are r = I 
and F = crT c 4 ff (see Sect. 4.1). 

It is worthy of note that boundary conditions at the 
disc surface have a very weak effect on midplane quantities 
as long as the disc is optically thick and has a large sur- 
face density, locally. Such an insensitivity is well known in 
stellar structure computations (Kippenhahn & Weigert, 
1990). Let us quote also that, contrary to a general be- 
lief, models that do not consider the atmosphere and use 
P = — as a boundary condition at z = h always vio- 
late the Eddington approximation since this relation does 
not guarantee that the atmosphere has the right optical 
thickness (i.e. I). Further, the approximation 



p(h) 



-dP 



P(oa) 9z 



-P 



Hz 



(18) 



becomes bad in regions where the surface density of the 
atmosphere is comparable to that of the disc, when the 
mean absorption within the atmosphere is low, or when 
the disc is not optically very thick. For these reasons, we 
sustain the idea that the use of the Eddington approxi- 
mation has sense only if the structure of the atmosphere 
is computed together with that of the disc interior. 



3. Ingredients and computational method 

3.1. Equation of state and opacities for a cosmic gas 

The physical problem depends on a certain amount of 
thcrmodynamical and radiative data. For the present ap- 
plication, the gas has cosmic abundances and is subject 



to atomic ionizations and molecular dissociations. Lo- 
cal Thermodynamic Equilibrium (LTE) is assumed (Hure 
et al., 1994b). The mean mass per particle /otih, coeffi- 
cients xt and Xpj adiabatic gradient V a d, heat capacity 
c p (needed to treat convective transport) and adiabatic 
exponent Ti have been computed accurately from chemi- 
cal abundances at thermal equilibrium (Hure, 1998). Grids 
of data with p and T as the inputs have been generated. 
Interpolations between mesh points are performed with 
bicubic splines. We have computed a high precision ana- 
lytical expression fitting p (and subsequently coefficients 
Xt and Xp) as functions of the temperature and density 
with an accuracy less than 4% relative to raw data in 
the case of a zero metallicity gas (see the Appendix, Sect. 
Bl). This fit can be used for a cosmic gas as well with- 
out producing important errors (since metals in a cosmic 
mix have very low abundances relative to hydrogen, with 
almost no influence on the EOS and related quantities). 

Opacity is probably the most important ingredient of 
the model and the choice for the grey absorption coef- 
ficient k is critical, specially near the surface (Mihalas, 
1978). Here, we take (Hameury et al., 1998) 



K = Oku + (1 — 9)kp 



(19) 



where kr and Kp are the Rosseland and Planck means 
respectively, and 9 is a function which varies continuously 
from to 1 as the optical depth increases from to infinity. 
A suitable choice is 



9m(r) 



1 



1 



(20) 



where the index m sets the stiffness of the transition from 
Kp to kb. (m = 1 is taken in the following applications). 
The ^-function is arbitrary and should strongly influence 
physical quantities at the disc surface. Also, we are aware 
that the best grey coefficient is probably neither given by 
Kr, nor by Kp, nor by Eq.([l9|). A flux weighted opacity 
mean could be better (Mihalas, 1978; Hubeny, 1990). 

Tables of Rosseland and Planck means have been taken 
from various sources (Pollack et al., 1994; Hure et al., 
1994b; Alexander & Fergusson, 1994; Henning & Stog- 
nienko, 1996, Seaton et al., 1994). But Planck absorption 
means published by the Opacity Project (OP) (Seaton et 
al., 1996) have not be selected due to apparently large er- 
rors of unknown origin (Alexander, 1998; Zeippen, 1999). 
Two grids of data have been generated after passing 
through filters (mainly running averages) to smooth out 
discontinuities. Finally, opacity grids cover the tempera- 
ture range 10 — 10 s K and extend from extremely high to 
extremely low densities where the ideal gas and LTE as- 
sumptions are expected to fail. An example of Rosseland 
and Planck means versus the temperature is displayed in 
Fig. (|l|) . As done for the equation of state and related quan- 
tities, opacity values between mesh points are interpolated 
using bicubic splines. Working with analytical opacities 
may present some advantages (Burgers & Lamers, 1989; 
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Fig. 1. Rosseland and Planck absorption coefficients ver- 
sus the temperature from various sources for hydrogen and 
metal mass fractions X — 0.70 and Z = 0.02 respectively 
and p = 10~ 15 x T 3 g.cm~ 3 (i.e. logi? = -3 in the OP 
format; see Seaton et al., 1994). Planck means have been 
shifted upwards for clarity. 



Collin & Dumont, 1990; Bell & Lin, 1994), but they can 
also produce numerical instabilities in the integration of 
the disc structure equations since they mostly consist in 
continuous but non derivable piecewise functions. 

3.2. Computational method 

The altitude of the top atmosphere being not known a 
priori, many numerical methods for solving Two Bound- 
ary Value Problems (TBVPs) are almost unusable in the 
actual situation. It is however possible to apply vari- 
able changes, namely V z — > V T in the atmosphere and 
V z — > V f in the disc interior (Cannizzo & Cameron, 1988) 
in order to recover a common TBVP with fixed boundaries 



(e-g- 



G [0, 1]), but this drastically increases vertical 



gradients near the surface and make the problem more un- 
stable from a numerical point of view. Some groups work 
with relaxation algorithms (Cannizzo, 1992; Milsom, Chen 
& Taam, 1994; D'alessio et al., 1998; Hameury et al., 1998) 
which require (good) vertical profiles as starting guesses 
and are known to be rapidly converging methods. How- 



ever, zones with steep gradients generally need a local 
mesh refinement. Other authors prefer straight integra- 
tors and algorithms based on shooting methods (Lin & 
Papaloizou, 1980; Meyer & Meyer-Hofmeister, 1982; Mi- 
neshige & Osaki, 1983; Smak, 1984; Mineshige, Tuchman 
& Wheeler, 1990; Rozariska et al., 1998; Papaloizou & 
Terquem, 1999) typical of Initial Values Problems (IVPs), 
with only a few quantities to guess but with possible trou- 
bles regarding the precision. 

In the present case, we have written a Fortran code 
named VS 3 KAD which performs the integration of the disc 
and atmosphere equations using single side shooting and 
numerical routines specially designed for stiff IVPs. Inte- 
gration steps are variable, internally to the routines and 
the accuracy can be as small as the machine precision. In 
its present version, the code works with the following vari- 
ables: lnT, hip, 



'"eff 



and r, but more clever choices 



are possible, specially to reduce artificially the stiffness of 
the equations. We have noticed that the computational 
time is considerably smaller than with Runge-Kutta algo- 
rithms, even with a variable step. Like in many stiff prob- 
lems, the numerical integration may turn out to be non 
conservative: integrating from the top atmosphere down 
to the midplane can give a solution which, in the surface 
neighborhood, can be very different than that obtained 
when integrating in the opposite direction (Mineshige, 
Tuchman & Wheeler, 1990; Press et al., 1992; De Kool 
& Wickramasinghe, 1999). This has been observed in the 
present case. It reminds the influence of both boundary 
conditions and underlying methodology. The best stabil- 
ity and reliability of the results are obtained by starting 
from the boundary layer where almost all the stiffness of 
the problem is concentrated. 

In practical, given a central mass M, accretion rate M, 
a-parameter, radius R and total energy deposition crT c 4 ff at 
the bottom atmosphere (see Sect. 4.1), the computation 
starts at an arbitrary altitude z = (the top atmo- 

sphere) where the surface density £( )(i?(°)) is guessed. 
The integration then proceeds downwards. Once the bot- 
tom atmosphere is reached, at an altitude z = h^°\ the 
equations for the disc interior are integrated down to the 
midplane where the net flux and surface density gener- 
ally differ from zero. By successive iterations on and 
'E(H^) — performed with a Newton- Raphson method — , 
quantities F«(0) and £^(0) can be driven to very small 
values. The problem has converged after n iterations, 
when simultaneously fW(0) w and E(")(0) w with 
the requested precisions 



fW(0) 



< € 



+ ' 



r 2E(")(0) y 



(21) 



(22) 
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and 



e H < 1 



< e 



+ ' 



(23) 



where e£, e+, e5, e_, e_ and are very small, positive 
values. Since our treatment of self-gravity is very approxi- 
mate, it is legitimate to allow a much lower (but sufficient) 
precision on the surface density (this accelerates conver- 
gence and lowers the computational time). On average, it 
takes a few milliseconds CPU on a single user personal 
computer and n < 25 with fo =e o = ( 6 o ) 2 = 10~ 10 
where = Sup(e F ,e^), e% = Sup(e s ,e^) and = 
Sup(e5,e^). This precision is sufficient in most cases. 
When self-gravity is neglected, Eq.(|J) becomes obsolete. 
The iteration is on the flux only and the code then runs 
much faster (by a factor ~ 5), with n < 6. 

The code has been tested and compared in many sit- 
uations (e.g. without self-gravity, without turbulent pres- 
sure, with and without external irradiation) with reference 
models (Rozariska, 1998; Hameury et al., 1998) and ap- 
pears to behave very well[jj It automatically scales to the 
central mass, and is able to model any kind of disc: AGN 
discs, CV discs, YSO discs and sub-nebulae as well. A sim- 
ulation of D/H enrichments in the PSN with this code is 
currently under way (Hersant, Gautier, Hure, 2000). 

4. A few applications of the model 

4-1- The hypothesis of a non irradiated disc 

Vertically averaged disc models show that the non self- 
gravitating and gas pressure dominated parts of ev-discs 
have an aspect ratio ^ ~ 0.001 — 0.1, depending on the 
central mass and accretion rate mainly. There is a slight 
flaring (i.e. jp-^ > 1) which eventually may vanish in the 
outermost regions (e.g. Ruden & Pollack, 1991; Bell et al., 
1997). The existence of a flaring inner disc is supported 
by observations: the spectrum of T-Tauri stars and AGN 
contains a prominent infrared component which is cur- 
rently interpreted as light re-processing in the superficial 
layers of the disc (Adams & Shu, 1986; Voit, 1991; Chi- 
ang & Goldreich, 1997). It is possible that self- irradiation 
also plays a role (Fukue, 1992). For circumstellar discs, 
the flaring angle of bare a-disc is too low to enhance the 
infrared spectral component to the observed level (Kenyon 
& Hartmann, 1987). And the heating by the central star 
should not produce a significant increase of the flaring (e.g. 
D'alessio et al., 1999). But the long wavelength emission 
can efficiently be boosted if the disc gets into a warped 
configuration (Terquem & Bertout, 1993; Miyoshi et al., 
1995; Bachev, 1999). 

Beyond a certain radius, self-gravity becomes impor- 
tant and reduces the disc thickness such that 



d In h 
d\nR 



< 



1 The author plans to make the executable available to the 
community. In the meanwhile, vertical structure computations 
may be performed on special request. 



(Sakimoto & Coroniti, 1981; Shore & White, 1982; Can- 
nizzo & Reiff, 1992; Liu, Xie & Ji, 1994; Hure et al., 1994a; 
Hure, 1998). As we shall see below, the same effect is ob- 
served from vertical structure computations. It means that 
outer parts should not receive directly light emitted at 
the center, contrary to the inner parts. Note that a disc 
is generally not isolated but embedded into a " warm" en- 
vironment which may in turn heat it up. Discs in YSO 
are surrounded by an envelope of gas and dust resulting 
from the cloud core collapse (Mundy, Looney & Welch, 
2000; Chick & Cassen, 1997; D'Alessio, Calvet & Hart- 
mann, 1997). In the case of AGN, clouds moving above the 
disc (BLR clouds) can scatter light back onto it (Shields 
1977; Collin-Souffrin, 1987; Osterbrock, 1993). It follows 
that, even outer regions which are not directly exposed 
to the luminous central source can be substantially irradi- 
ated. The disc response to irradiation is a complex prob- 
lem to solve self-consistently since it depends on many 
parameters (size and location of the ionizing source(s), 
shape of the ionizing spectrum, intensity of irradiation, 
disc flaring angle, surface albedo, disc optical thickness, 
gas metalicity, etc.) which are neither well known, nor 
well constrained by current observations. For discs having 
a large total surface density, mostly the superficial layers 
are heated up and ionized, and the midplane is essentially 
not affected (Sincell & Krolik, 1997; Collin & Hure, 1999; 
Nayakshin, Kazanas & Kallman, 1999; Igea & Glassgold, 
1999). Deep structural changes (for instance, a tempera- 
ture inversion or an iso-thermalization) may however oc- 
cur in some cases (D'alessio et al. 1998). Although discs 
should experience some external heating, even in their out- 
ermost (self-gravitating) parts, we do not consider irradi- 
ation in this paper, for simplicity. So, considering internal 
viscous heating only, the disc effective temperature, far 
from the central object, is given by 



8 7TCT 



-n 2 M 



(24) 



where M is the mass accretion rate. 

In the following, we apply the code to the computation 
of the vertical structure of two different systems subject 
to vertical self-gravity: circumstellar discs and AGN discs. 
We restrict to central masses of 1 M Q and 10 s M Q respec- 
tively, and we show on a few examples how self-gravity 
and turbulent pressure which are the main specificities of 
this model affect the disc structure. Computations are per- 
formed discarding any kind of instabilities the disc could 
undergo. These are stopped when the temperature at the 

top atmosphere (i.e. T(H) — (j) 1 ^ 4 T c g) attains 10 K, for 
both physical and practical reasons. 

4-2. Effect of turbulent pressure on the vertical 

temperature and density stratification: an example 

To see clearly the influence of turbulent pressure (see 
Eq.([l0|)), we have chosen a simulation with a = 1 which 
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13 
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Fig. 2. Effect of turbulent pressure on the distribution 
of the density (top) and temperature (bottom) with the 
altitude z in an AGN disc with M = 10 8 M , M = 1CT 1 
M Q /yr and a = 1 at 500 and 1000 i?s from the center (Rs 
is the Schwarzchild radius of the black hole) . The position 
of the bottom atmosphere is marked by a circle. 



corresponds to the very upper limit regarding supersonic 
turbulence. In addition, convection and self-gravitation 
are turned off. Under these circumstances, we have com- 
puted the vertical structure with viscosity law v\ for an 
disc around a massive black hole with M = 10 s Mq, at 
two different radii: R = 500 R s (R s = 2GM/c 2 is the 
Schwarzchild radius of the black hole) and R = 1000 Rs- 
The accretion rate is M = 10 _1 M Q /yr. The temperature 
and density profiles so obtained are plotted in Fig.(|J). 
We see that turbulent pressure makes the disc thicker, as 
expected intuitively. Actually, according to Eq.(|IT|), the 
density gradient is lowered with turbulent pressure and 
the transition from the ambient medium to the disc in- 
terior is softer. The flux gradient depending on the den- 
sity to a positive power — in the absence of convection at 
least — , a larger integration range is needed to obtain a 
zero net flux at the midplane, hence a thicker disc. The 
effect is specially visible at R = 1000 Rs, with a 41% disc 
thickening (measured on h; we have 33% on H). The disc 
flaring therefore increases by the same relative quantity. 
Since we have not explored the whole parameter space, 
it is probably easy to find cases where the effect is more 



important. The total surface density is almost unchanged 
(a 2% increase only). At 500 Rs, the disc thickening is less 
important, about 15%. Note in this example the presence 
of a density inversion which occurs well below the base of 
the atmosphere. It tends to be washed out by turbulent 
pressure, as argued in Sect. 2. 3. 

The way turbulent pressure acts on both the tempera- 
ture and density at the midplane temperature is however 
less straightforward and depends intimately on the verti- 
cal stratification. For a given effective temperature, a disc 
with a larger surface density would theoretically have a 
hotter core. But this happens only in the calculation at 
500 Rs, and the effect is very minor. 

We have checked that the conclusions derived here- 
above are qualitatively unchanged using 1/2 (see Sect. 4.5 
for another effect of turbulent pressure) . We will not dis- 
cuss the case of circumstellar discs because values of the a- 
parameter currently adopted for these systems are rather 
very low (typically a ~ 10~ 3 — 10~ 2 ; D'Alessio, Calvet & 
Hartmann, 1997) and so, no effect of turbulent pressure is 
indeed observable. This has been verified. 

4-3. Effect of the viscosity law: v\ versus V2 

In this paragraph, we show a few differences between v\ 
and vi. Convection is taken into account and turbulent 
pressure is left aside. The midplane temperature, total 
surface density, disc thickness and disc mass are plot- 
ted versus radius in Fig.(|) for M = 1 M , M = 10~ 7 
MQ/yr and a — 10~ 3 . These parameter values are typ- 
ical of discs around T-Tauri stars (D'Alessio, Calvet & 
Hartmann, 1997). The disc mass M disc is 



M disc (R) = 2tt / R'T, t (R')dR' ^ nZ t R 2 



(25) 



where i?j n is the inner edge of the disc (unimportant as 
long as R ^> R[ n )- Results computed without self-gravity 
are also shown in comparison. We see that, in the non self- 
gravitating region, differences on the geometrical thick- 
ness and midplane temperature are small, the disc be- 
ing slightly thicker and hotter with v\ than with V2- The 
most important effect is on the surface density and con- 
sequently on the disc mass: the ^i-prescription leads to 
slightly more massive disc than the ^-prescription (by a 
factor ~ 2 here). 

We see that the disc is definitely affected by self-gravity 
from approximately 7 AU whatever the viscosity law (in 
fact, a little bit less with v\ which agrees with the fact 
that that prescription leads to more massive discs, as long 
as M dlsc < M). Let us remind that the importance of 
self-gravity can be measured by the quantity (see Eq.(||)) 

47rGE(i?,z) 



C(R,z) = 



n 2 z 



(26) 



which remains finite at z — as a Taylor expansion shows. 
It follows that a " natural" definition for the limit between 
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the standard disc and the self-gravitating disc can be £ = 1 
at the midplane (if £ > 1 then the disc contribution to 
gravity exceeds that due to the central object, vertically). 

On the example, the disc thickness goes through a 
maximum near 10 AU for both viscosity laws. Beyond this 
radius, the disc becomes thinner and thinner. Regions lo- 
cated farther away can therefore not intercept photons 
emitted at the center. Interestingly, the temperature and 
the surface density show very different behaviors. With 
law v\ , the temperature decreases monotonically as R in- 
creases, almost as in the absence of self-gravity. So does 
the total surface density. The midplane density, not shown 
on graphs, gently increases with the radius. Note that St 
(and Af dlsc ) are surprisingly not affected by self-gravity. 
The origin of this insensitivity has not been identified at 
the time being. Conversely, with law v^, the surface den- 
sity and the density violently increase with the radius. 
The disc thickness falls in also very rapidly. The temper- 
ature reaches a plateau. In fact, this somewhat "singular" 
behavior is predicted by the vertically averaged model 
(Shlosman & Belgelman, 1987; Hure et al. 1994a; Hure, 
1998; see also Duschl, Biermann & Strittmatter, 2000). In 
particular, the temperature T sg at the plateau is mainly 
fixed by the ratio M/a, namely 



T — 

J-sg — 



4G 2 M 2 



1/3 



9 a 2 



24100 n 



M 



2/3 



1 M /yr 



-2/3 



K 



(27) 



For M/a = 10" 4 M /yr as in Fig.(g), Eq.@ yields 
T sg ~ 120 K (assuming [i — 2.33). Although this value 
is derived from the one zone model, it is in good agree- 
ment with the vertical structure computation which gives 
a plateau at about 70 K. Note that the drastic increase 
of the surface density produces an important increase of 
the disc mass which even exceeds the value obtained with 
law v x . When M disc > M (this occurs at 23 - 28 AU 
on the example), the disc becomes very self-gravitating, 
probably gravitationally unstable (Goldreich & Lynden- 
Bell, 1965; Shu et al., 1990) and should therefore not be 
well described with current viscosities (e.g. Lin & Pringle, 
1987) nor by steady state solutions (Goldreich & Lynden- 
Bcll, 1965). One reaches here the limit of the model. It 
is however interesting to note that, asymptotically, there 
is either a hot very dense solution or a cold more diffuse 
solution (see Paczynski, 1978). 

The same quantities obtained under the same condi- 
tions but for an AGN disc with M = 10 8 M , M = 10~ 2 
M /yr and a = 0.1 are displayed in Fig.(^). Globally, we 
find similar trends. The singular behavior observed in the 
previous example with v<i inside the self-gravitating regime 
seems even stronger. Self-gravitation becomes important 
from about 500 Rg. The midplane temperature stabilizes 
radially at ~ 3900 K. This value again compares quite well 



a 3 - 



o 2 - 



14 





\\ 


- 


\\ 




\\ 












v\ 




\\ 




\\ 




\\ 




\\ 




A 




\ 


1 1 



13 




— v,, no self-gravity 

— with self-gravity 

— v 2 , no self-gravity 

— with self-gravity 
#C(R,0)=1 






log R A( 



Fig. 3. Midplane temperature, total surface density, geo- 
metrical thickness and disc mass versus the radius (in AU) 
computed with v\ (solid lines) and V2 (dashed lines) for a 
disc with M — 10~ 7 M /yr and a = 10~ 3 around a one 
solar mass star. 



with Eq.@ which gives 6600 K (assuming fi = 1.27). The 
presence of extremely steep surface density and density ra- 
dial gradients probably means that the assumption made 
on the disc potential (see Sect. 2.1) is no longer valid. This 
would be worthwhile to check. The physical picture there 
is then that of a disc surrounded by a very dense ring 
(E increases by more than one order of magnitude over 
a relative length ~ 15%) which contains almost all 
the disc mass and angular momentum. When M dlsc > M 
(beyond a few 10 3 Rs for law z^), effects of self-gravity are 
expected to be global, with a change in the rotation law. 

All these results require a few comments. First, it is 
important to note that the effect of self-gravity is slightly 
underestimated when it is neglected in the computations. 
This is specially true with law V2- Besides, self-gravity 
becomes important before £ reaches unity, like in the one 
zone model (Hure, 1998). At the inner edge of the self- 
gravitating disc (that is at C = 1 following our definition), 
the disc mass is much less than the central mass: we find 



M dlsc 
M 



io- 



10 1 (depending on the viscosity law) in 



the YSO case, and 



M dl3c 
M 



1Q- 3 in the AGN case. This 



suggests that low mass discs must not be automatically 
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Fig. 4. Same legend as for Fig.(||) but for an AGN disc 
with M = 10 8 M Q , M = 10~ 2 M /yr and a = 0.1 (R s is 
the Schwarzschild radius of the black hole) . 

classified as non self-gravitating discs as often asserted. In 
particular, it is not excluded that T-Tauri discs we observe 
be subject to vertical self-gravity, despite their relatively 
low mass (Beckwith et al., 1990). 

It has been stated in Sect. 2.4 that the relation be- 
tween v\ and V2 is not trivial. In the non self-gravitating 
parts however, midplane temperatures T(0) show a strik- 
ing quasi-parallel variation with the radius in logarithmic 
scales (this is also true for quantities E t and h), which 
could be attributed to the fact that v\ is indeed propor- 
tional to V2 ( or one goes from one solution to the other 
by a change of a; see power law solutions for a-discs). 
To demonstrate that these two viscosity prescriptions are 
intrinsically different (as suggested by what happens in 
the self- gravitating regions), we have computed the ra- 
tio j£ = |f2A p ri/c s with v\, for the examples discussed 
above. We have chosen two radii, the one lies in the classi- 
cal disc part, the other is in the self-gravitating part. The 
results are shown in Fig. (^) . 

4-4- Example of internal structure: 2D-density maps 

We show in Fig. (j^) the density field within a circumstel- 
lar disc for the same parameter values as in Fig.(^). This 




z/h 



Fig. 5. Ratio v\/v2 versus z/h using V\ for an YSO disc 
(same parameter values as for Fig.([|)) at 1 AU (top left) 
and 10 AU (bottom left), and for an AGN disc (same pa- 
rameter values as for Fig.(Q)) at 100 Rs (top right) and 
500 Rs (bottom right). 



simulation has been carried out using V\ and v 2 as well, 
including convection, self-gravity and turbulent pressure. 
We also indicate iso-contours of the temperature, the limit 
between gas pressure dominated zones and radiative pres- 
sure dominated ones, and lines where C,{R, z) = jk, 1 and 
10. We notice that the density and the temperature show 
strong variations in between the midplane to the photo- 
sphere's base and we are far from a vertically isothermal, 
homogenous disc as often considered. Note, specially for 
viscosity law u\, the re-increase of the midplane density 
in the radial direction (another density inversion) as self- 
gravity gains in importance. It is remarkable that self- 
gravity does not appear suddenly (i.e. on a small radial 
range) but installs gently and steadily over a very ex- 
tended domain. For instance, with law i>2, it contibutes 
by 10% in the hydrostatic equilibrium near 2 — 4 AU and 
by 50% at about 7 AU, and all disc quantities are signifi- 
cantly modified over this domain. We see clearly that the 
disc density is slightly smaller with v 2 than with v\ (same 
color code for both maps). 

We notice that the whole disc is optically very thick, 
due to its large surface density. It is then likely that a 
mean external heating should not change noticeably the 
temperature and density at the midplane. From this point 
of view, the neglect of irradiation is entirely justified. It is 
interesting to see that, as long as an a-model can be ap- 
plied to the Solar Nebula (e.g. Ruden & Pollack, 1991; Pa- 
paloizou & Terquem, 1999; Drouart et al, 1999), our giant 
planets (displayed on graphs) lie inside the self-gravitating 
part of the disc (this is true whatever the accretion rate, 
but depends on the a-parameter; see the next Section). 
As already noticed (Ruden & Pollack, 1991), this coinci- 
dence is somewhat striking and we do not know to what 
extent self-gravity could have played a role in the process 
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Fig. 6. Density field within a disc surrounding a one solar mass central star, computed with v\ (left) and with 
(right). The accretion rate and viscosity parameter are M = 10 -7 M Q /yr and a = 10~ 3 respectively. The color code 
refers to logp (with p in g/cm 3 ) and is the same for both maps. Iso- values of logT are given in dotted lines. Also 
indicated in bold lines are contours for ( = 0.1, 1 and 10. The limit between the radiative pressure dominated zone 
(RPDZ) and the gas pressure dominated zone (GPDZ) is marked by a thin plain line. 
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of the formation of outer planet and other objects (Lis- 
sauer, 1993). 

Similar 2D-maps are displayed in Fig.(Q) for an AGN 
disc with M = 10 8 M 0! M = 10~ 2 M Q /yr and a = 10" 1 , 
using the same color table. Note the existence of an in- 
ner radiative pressure dominated region at R < 100 i?s, 
common for AGN discs as well as the great steepness of 
density gradients near the surface which may cause nu- 
merical difficulties. Also visible is the density inversion at 
R ~ 160 — 200 Rs which has less amplitude and is less 
extended radially in the presence of convection and tur- 
bulent pressure. 

4-5. A criterion to check the importance of vertical 
self-gravity 

We see on both Figs.(§) and (g) that ((R,oo) < C(R,0) 
and iso- values of £(i?, z) form very curved lines ( "diabolo"- 
shape surfaces in 3D). This is important if one wishes 
to estimate correctly the position of the self-gravitating 
regime. For instance, in the case of the law v\ discussed in 
Fig.(||), £ = 1 in the midplane at 7 AU (we have 450 Rs 
in the AGN case depicted in Fig.(0)) whereas this occurs 
in the disc atmosphere at about twice the distance. 

For many purposes, it is interesting to know the posi- 
tion of the (vertically) self-gravitating disc. That is why 
we have performed a systematic computation to find the 
radius R sg satisfying ((R sg ,0) = 1 as a function of the ac- 
cretion rate, for 4 values of the a-parameter in the range 
10~ 3 - 1, and for M = 1 M and M = 10 s M Q . These 
parameter values should cover a wide variety of circum- 
stellar discs and AGN discs. The calculations have been 
performed with vi including convection, self-gravitation, 
with and without turbulent pressure (the use of v\ would 
systematically give a slightly lower value of R sg ). The re- 
sults are plotted in Fig.(||) for parameter values typical 
of circumstellar discs. We see that the location of self- 
gravitating regime is very sensitive to the viscosity pa- 
rameter. The larger the value of a, the further away the 
self-gravitating region. Regarding the sensitivity to the ac- 
cretion rate, there are three different trends. For a < 10~ 2 , 
the lower M, the larger R sg . For a ~ 0.1, the dependence 
is rather weak: R sg « 27 ± 3 AU. For higher values of 
the viscosity parameter (generally not appropriate to fit 
properties of observed discs), R sg is an increasing function 
of M. Note also that turbulent pressure, when important, 
slightly increases R sg (~ 12% in the example for a = 1 
and the highest accretion rate). 

Results for AGN discs are shown in Fig.(^). As above, 
large viscosity parameters imply wider non self-gravitating 
inner disc (because discs are less dense with large values 
of the a-parameter), and turbulent pressure pushes fur- 
ther away the self-gravitating regime. For a < 10 -2 , we 
have R sg ~ 5 — 20 x 10 15 cm. The sensitivity to the accre- 
tion rate remains weak. For higher values of the viscosity 
parameter and for very low to moderate accretion rates 
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Fig. 8. Inner edge of the self-gravity dominated disc R sg 
versus the accretion rate for a = 10 -3 , 10~ 2 , 0.1 and 1 
with (bold lines) and without turbulent pressure (thin 
lines). Lines of constant effective temperatures are in 
dashed lines. 
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Fig. 9. Same legend as for Fig.(|) but for M = 10 8 M Q . 

(M < IQ- 2 - 10- 1 M Q /yr, depending on a), R sg glob- 
ally decreases as M increases. For higher accretion rates, 
i? sg increases as M increases and self-gravity is important 
within the radiative pressure dominated region. 

It is important to notify that the "Toomre's criterion" 
(Qt = c *q^ < 1) is an helpful tool to trace regions 
where some gravitational instabilities occur in a stellar 
disc (Toomre, 1965), but it differs significantly from the 
criterion derived by Goldreich & Lynden-Bell (1965) for 
gaseous discs (Qglb = ^ 1)- Besides, Qt is com- 
monly written in different forms which are not strictly 
equivalent (e.g., Ruden & Pollack, 1991; Sincell & Krolik, 
1997; D'alessio, Calvet & Hartmann, 1997, Papaloizou & 
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Fig. 10. Comparison between quantity £ at three differ- 
ent altitudes around the transition C = 1 (computed with 
self-gravity) and the Toomre's parameter Qt computed 
(without self-gravity) for three different definitions, for an 
AGN disc (top; same parameter values as in Fig.®) and 
an YSO disc (bottom; same parameter values as in Fig.(^)). 

Terquem, 1999), specially in a 2D- model. Althought these 
Q-parameters are undisputably related to the ^-parameter 
in some ways (for instance, Qt X C(R> °°) = ^qh > 1 
and Qglb x ({R, 0) ~ 16), they must be used with care if 
one wishes to check the importance of self-gravity. They 
may lead to noticeable over-estimates of R sg , specially if, 
as it is always the case, the Qx-parameter is computed 
from models that discards self-gravity. This is illustrated 
in Fig. (jT^) where we have plotted values of C a t three 
key altitudes and Qt versus the radius. It seems therefore 
preferable to make the check on ((R, 0) rather than on any 
other quantities, otherwise the importance of self-gravity 
is under-estimated. 

5. Conclusion 

In this article, we have presented the equations account- 
ing for the vertical structure of steady state keplerian 
accretion discs, including simultaneously turbulent pres- 
sure, convective transport in the framework of the Mixing 
Length Theory, and the disc self-gravity within the infi- 
nite slab approximation. The emphasis has been placed on 
outer regions where vertical self-gravity becomes increas- 
ingly influencal. We have shown that turbulent pressure 
makes the disc thicker, but it can be neglected if a < 0.1. 
Also, self-gravity may be important even if the disc mass 
is very low compared to the central object, contrary to 
usual assertions. Further, the transition from the classical 
disc to the self-gravitating disc takes place gradually and 
so, concerns a large radial domain. 



Another important conclusion of this work is that the 
model of mechanical energy deposition towards the ver- 
tical direction is crucial. This is not surprising since vis- 
cosity is the source of heating and accretion. Two differ- 
ent behaviors are predicted inside the self-gravitating (and 
gravitationally unstable) region. On the one hand, the otP- 
formalism yields a solution which is asymptotically cold 
and diffuse, whereas the local version of the standard pre- 
scription gives a more " singular" solution, asymptotically 
hot and dense, which is also predicted by the vertically av- 
eraged model. We can imagine that a suitable choice for 
the function v(z) should provide any intermediate solu- 
tion, meaning that this model for outer regions has almost 
no predictive power. Although the present stationary ke- 
plerian disc model is probably irrelevant to describe the 
disc structure in its strongly self-gravitating parts, the ex- 
istence of two asymptotically distinct solutions might in- 
dicate that the outer disc can get into, either a flat and 
cold configuration where fragmentation and formation of 
indiviual clouds (Kumar, 1999) and compact objects like 
planets and stars (Collin & Zahn, 1999) could occur, or 
into a thick diffuse configuration (a torus) (Paczynski, 
1978). Time dependent simulations should shed light on 
this question. 

This study confirms that discs around super-massive 
black holes are self-gravitating close to the center, be- 
yond a few hundreds Schwarzchild radii, depending on 
the central mass accretion rate and a-parameter. We have 
found that the disc surface density remains high when self- 
gravity dominates and that the disc mass definitely rises 
outwards. Let us remind that a current problem with the 
fueling of active nuclei is that the accretion time scale 
M C /M is usually much too long at large radii, that is 
why other more efficient mechanisms are invoked (Frank, 
1990). However, given the sensitivity of the disc struc- 
ture to viscosity, it is not excluded that depth dependent 
viscosity laws other than considered here would lead to 
smaller surface density distributions, and consequently to 
shorter accretion time scales. It is therefore important to 
test other models for the energy deposition along the z- 
axis as well as other kinds of viscosity prescriptions (Lin 
& Pringle, 1987; Cannizzo & Cameron, 1988; Richard & 
Zahn, 1999; see Hure & Richard, 2000). 

With accretion rates and values of the viscosity param- 
eter usually considered to model discs in T-Tauri stars, we 
conclude that our giant planets (if not Jupiter, the other 
ones) were probably formed within the self-gravitating 
part of the Solar Nebula (Drouart et al. 1999). This is 
a fortiori true if these objects experienced any inward mi- 
gration. More generally, as long as the a-theory may be 
applied to describe the inner parts of circumstellar discs 
around forming stars, self-gravity is expected to play a role 
at about ten AU from the center. This is not uncompatible 
with the discs we observed (Beckwith et al., 1990). 
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Appendix A: Treatment for convection in the 
presence of self-gravitation 

Convection is treated with in the framework of the Mixing 
Length Theory (MLT) (Cox & Giuili, 1968). It is modified 
in order to include turbulent pressure in the determina- 
tion of the pressure scale height, and self-gravititation. 
Following the original version of the MLT, we neglect the 
variation of gravity with the altitude z. It is likely that 
turbulent pressure modify significantly the velocity of ris- 
ing elements within convectively unstable zones, specially 
if the a-parameter is close to unity, but this effect is not 
considered here. The total pressure height scale writes 

p d\n{P+p t ) y ' 

where P is the total (gas plus radiation) pressure and p t 
is the turbulent pressure. This expression being singular 
at the midplane (X p — > oo) (see Eq.(|])), the pressure scale 
height is usually limited to the disc thickness 

P + Pt 



A p = Inf 



(A.2) 



p(f7 2 z + 47rGS)_ 
where p is the gas density, f£ is the rotation frequency, E 
is the surface density and h is the altitude of the bottom 
atmosphere. 

According to Eq.(Q), the temperature gradient is 

iJ sr = -f < A ' 3 > 

dZ An 



2 r p depends on the adia- 
dinP7ad ( com P u ted from the EOS; 



where the actual gradient V 
batic gradient V a d = 1 " ' 

see the next Section) with respect to the fictitious radia- 
tive gradient V r defined as 

3 pnFXp 



Vr 



(A.4) 



16 crT 4 

where F is the flux to be transported upwards and k 
is a grey absorption coefficient which must approach the 
Rosseland mean kr. at great optical depth (see Sect. 3.1). 

When V r < V ac i (i.e. in convectively stable zones), 
heat is transported through radiation only and we have 

V = V r - Conversely, when V r > V a d, V is computed 
from 

V = (1 - z 3 )V r + x 3 V ad (A.5) 
where x is the real root of the third degree equation 
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where ckmlt is the mixing length parameter («mlt = 1-5 
here) and c p is the constant pressure specific heat capacity. 
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Fig. A.l. Thermodynamic coefficients /i, xpi XT, c p , T\ 
and V a< a computed at thermal equilibrium for a hydrogen- 
helium mixture with elementary abundances H:He=l:^j 
for log/? = —12 (plain lines), log/9 = —10 (dotted lines) 
and log/? = —8 (dashed lines). Raw data are in bold. 
For n, Xp and xt, the difference between raw data and 
the high precision fitting formula given in the Appendix 
(Sect. Bl) is not visible. For all coefficients, thin lines are 
obtained from low precision approximate expressions (see 
the Appendix, Sect. B2). 



Appendix B: Notes on the EOS and related 
quantities 

B.l. High precision fitting formula for // as functions of 
the temperature and density — x~ coe ffi c i en t s 

Within the assumption of LTE, the pressure of a mixture 
of radiation and ideal gas undergoing atomic ionization 
molecular dissociation is given by Eq.(j^). In particular, 
the mean mass per particle fi, in units of the proton mass, 
is defined as 



J2i m iPi 



(B.l) 



where and Pi are respectively individual mass and par- 
tial pressure of the chemical compounds, and the sum 



Table B.l. Coefficients and functions required to com- 
pute n(p, T) according to Eq.(B.2). Note that S1+...+S4, = 
0.618 (fully ionized helium), and So — (Si + •■■ + ^4) = 2.373 
(fully molecular gas). 
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extends over the total number of compounds, including 
electrons. We have computed equilibrium abundances for 
a mixture of hydrogen and helium (H:He=l: jq) only as 
function of p and T (see Hure, 1998). Data relative to the 
EOS so obtained are plotted versus the temperature in 
Fig. (A.l) for three values of the gas density. In particular, 
the resulting function /i(p, T) has then been fitted with 
great accuracy by the following expression 



fit 



with 



So 



i=l,4 



Si tanh 



iog(r) - 



(B.2) 



(B.3) 



where coefficients Si and and Aj are listed in Tab. flB.l 
Hyperbolic functions account successively for the transi- 
tions H2/HI, HI/HII, He/Hel and Hel/Hell. Relative er- 
rors on raw data —fr = M ,7 M are displayed in Fig. (B.l) 



and never exceed 4 % over the whole domain of temper- 
ature and density. Besides, the standard deviations with 
respect to raw data are less than 1.3 %. We have tried to fit 
the residuals with a few Gaussians, but no major improve- 
ment has been obtained. Note that, by construction, the 
fitting formula shows no singular behavior at high of low 
densities and temperature. The fit can be used for a cos- 
mic gas without producing large errors (the effect of heavy 
elements on the EOS is very weak as long as Z <C X). 

The temperature and density exponents of the total 
pressure are repectively given by 
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Fig. B.l. Relative error on the mean mass per particle 
between raw data and fitted values, versus logT (upper 
panel) and versus logp (lower panel). The standard devi- 
ation is plotted in bold line. 



and 

X P = 



dlnP 
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where is the ratio of gas pressure to total (gas plus 
radiation) pressure. It follows from Eq. 
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and the precision with respect to raw data is also of the 
order of a few percents. 



where U is the specific internal energy of gas and radia- 
tion, 



c XtP 

P XpPT^ad 

and 
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Since disc models are still very uncertrain, one may 
use, as a first approximation, low precision expressions for 
all coefficients related to the EOS discarding changes in 
the internal energy due to ionizations and dissociations, 
for instance 



Xp~P 
XT ~ 4 - 3/3 
and then 
V ad ~ 
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3/3 
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B.2. Low precision expressions for coefiicients \t, Xp> 
V ad , r x and c p 

The adiabatic gradient V a d, heat capacity at constant 
pressure c p and adiabatic exponent Ti are respectively 
defined by 



